BARRETO et al. 2022: Habitat loss and biodiversity gain (BIOCON-S-22-00746)

Model formulas

For each of our response variables we run three candidate models for each scale:

Model 0 <- ~1, absence of effect

Model 1 <- percentage of forest cover at 3 and 5km radius

Model1.2 <- non-linear forest cover at 3 and 5km radius (quadratic model)

Loading data and gathering into one datasets calculated in script ‘data_wrang_GIT.Rmd’, content: 1) “beetles”: raw data for species occurrence per landscape/site; 2) “hab.data”: data containing classes of habitat association; 3) “env.data”: environmental variables, forest cover percent at 3 and 5km radius buffer 4) “div_hab”: diversity (alpha, beta and gamma) between classification habitat association 5) “rc.df”: Beta RC calculations 6) “ab_hab”: mean and median abundance between groups of habitat association 7) “st.hab_ab”: abundance data between classification habitat association per site

Landscape gamma_FS gamma_NFS alpha_FS alpha_NFS betad_FS betad_NFS ab.land_FS ab.land_NFS brc.fsNA brc.nfsNA brc.abund.fsNA brc.abund.nfsNA fc_3km fc_5km
148 11 9 3.75 2.25 7.25 6.75 66 40 -0.2917857 0.1800000 0.0007150 0.4310978 38.52523 43.86337
215 19 7 5.62 2.00 13.38 5.00 174 51 0.1050000 -0.3104762 0.2842485 0.3969207 30.15485 25.23160
263 16 13 6.00 4.50 10.00 8.50 393 214 -0.4914286 -0.0328571 0.6290100 0.3908194 10.10641 13.41180
266 11 8 3.25 1.38 7.75 6.62 70 12 -0.2290476 0.1309524 0.1069164 0.1522952 30.87525 31.16090
282 17 11 6.50 4.12 10.50 6.88 219 212 -0.3271429 -0.1109524 0.2140712 0.7270127 15.26270 16.35308
291 9 3 3.71 1.71 5.29 1.29 150 49 -0.3876190 -0.5861905 0.0977168 -0.4293818 45.74294 45.13213
317 12 4 5.00 1.38 7.00 2.62 152 56 -0.3832143 -0.5346429 0.1839339 -0.4355784 48.75805 47.85138
329 21 14 8.62 5.75 12.38 8.25 369 262 -0.0821429 0.2232143 0.9174174 0.4880952 14.73791 13.62588
333 14 15 4.50 4.62 9.50 10.38 158 219 -0.1314286 0.1921429 0.6792864 0.6124339 13.71403 14.50960
335 15 8 5.75 2.38 9.25 5.62 140 113 -0.2553571 -0.3357143 0.3477763 -0.9467563 23.87913 19.21670
359 16 17 6.00 5.50 10.00 11.50 149 133 -0.2253571 0.3067857 0.7706278 0.8152081 18.10208 17.62431
399 20 14 7.88 5.12 12.12 8.88 432 169 -0.1242857 0.2692857 0.9614376 0.7002717 29.73616 24.97851

RESULTS

We found 28 species to be forest specialists, those that come from the Atlantic Forest domain and exclusively collected in forested habitats; while 23 were non-forest specialists, species from a broader domain and/or less restrictive to forest habitats (hereafter FS and NFS, respectively).

FOREST SPECIALISTS

Gamma diversity

Modelling

g0 <- glm(gamma_FS ~ 1, data=data_div, family= Gamma (link = "identity"))
g5k <- glm(gamma_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
g5kq <- glm(gamma_FS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
g3k <- glm(gamma_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
g3kq <- glm(gamma_FS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))

## AIC-based model selection
(g.fs.aic <- AICctab(g0, g5k, g5kq, g3k, g3kq,
                      base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## g5k  -27.9   64.8   4.7     0.0 3  0.731 
## g3k  -29.6   68.2   3.0     3.4 3  0.136 
## g5kq -27.8   69.4   4.8     4.6 4  0.074 
## g0   -32.6   70.6   0.0     5.8 2  0.041 
## g3kq -29.2   72.2   3.4     7.3 4  0.019

Summary top model

summary(g5k)
## 
## Call:
## glm(formula = gamma_FS ~ fc_5km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.23181  -0.12743  -0.03948   0.15159   0.27856  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 20.68384    1.96964   10.50 1.01e-06 ***
## fc_5km      -0.21451    0.05782   -3.71  0.00404 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.0352002)
## 
##     Null deviance: 0.74994  on 11  degrees of freedom
## Residual deviance: 0.34328  on 10  degrees of freedom
## AIC: 61.839
## 
## Number of Fisher Scoring iterations: 3
g.fs <- g5k# rename top model to save and plot below

Inspect residuals:

sim <- simulateResiduals(g5k)
plot(sim)

testDispersion(sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.1487, p-value = 0.632
## alternative hypothesis: two.sided
hnp::hnp(g5k)
## Gamma model

boot::glm.diag.plots(g5k) ## for glm poisson or neg binom

Prediction:

Alpha diversity

Modelling

a0 <- glm(alpha_FS ~ 1, data=data_div, family= Gamma (link = "identity"))
a5k <- glm(alpha_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
a5kq <- glm(alpha_FS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
a3k <- glm(alpha_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
a3kq <- glm(alpha_FS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))


## AIC-based model selection
(a.fs.aic <- AICctab(a0, a5k, a5kq,
                     a3k, a3kq, 
                     base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## a5k  -19.3   47.6   2.8     0.0 3  0.523 
## a3k  -20.3   49.5   1.9     1.9 3  0.199 
## a0   -22.1   49.6   0.0     2.0 2  0.195 
## a5kq -19.0   51.8   3.1     4.2 4  0.064 
## a3kq -20.3   54.2   1.9     6.6 4  0.019

Summary top model

summary(a5k)
## 
## Call:
## glm(formula = alpha_FS ~ fc_5km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.43224  -0.11774  -0.03674   0.08987   0.35754  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.42155    0.94502   7.853 1.38e-05 ***
## fc_5km      -0.07200    0.02815  -2.558   0.0285 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.06034959)
## 
##     Null deviance: 0.97549  on 11  degrees of freedom
## Residual deviance: 0.61269  on 10  degrees of freedom
## AIC: 44.608
## 
## Number of Fisher Scoring iterations: 4
a.fs <- a5k# rename top model to save and plot below

Inspect residuals:

sim <- simulateResiduals(a5k)
plot(sim)

testDispersion(sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.1043, p-value = 0.656
## alternative hypothesis: two.sided
hnp::hnp(a5k)
## Gamma model

boot::glm.diag.plots(a5k)

Prediction:

Beta diversity additive

Modelling

### LM Beta for regular model selection
b0 <- glm(betad_FS~ 1, data=data_div, family= Gamma (link = "identity"))
b5k <- glm(betad_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
b5kq <- glm(betad_FS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
b3k <- glm(betad_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
b3kq <- glm(betad_FS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))

## AIC-based model selection
(b.fs.aic <- AICctab(b0, b5k, b5kq, 
                     b3k, b3kq, 
                     base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## b5k  -22.1   53.3   5.1     0.0 3  0.690 
## b5kq -21.4   56.6   5.8     3.3 4  0.133 
## b3k  -23.9   56.8   3.3     3.6 3  0.116 
## b3kq -22.7   59.2   4.5     5.9 4  0.036 
## b0   -27.2   59.8   0.0     6.6 2  0.026

Summary top model

summary(b5k)
## 
## Call:
## glm(formula = betad_FS ~ fc_5km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.24293  -0.12664  -0.05680   0.08772   0.34333  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 13.29031    1.24634  10.663  8.8e-07 ***
## fc_5km      -0.14349    0.03623  -3.961  0.00268 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.03508083)
## 
##     Null deviance: 0.76409  on 11  degrees of freedom
## Residual deviance: 0.32769  on 10  degrees of freedom
## AIC: 50.257
## 
## Number of Fisher Scoring iterations: 4
b.fs <- b5k # rename top model to save and plot below

Inspect residuals:

sim <- simulateResiduals(b5k)
plot(sim)

testDispersion(sim) # ok

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.1585, p-value = 0.608
## alternative hypothesis: two.sided
hnp::hnp(b5k) ## for glm Gamma or neg binom (maybe also others)
## Gamma model

boot::glm.diag.plots(b5k) ## for glm Gamma or neg binom

Prediction

bRaup-Crick (presence/absence-based)

Modelling

### LM gamma for regular model selection
rc0 <- lm(brc.fsNA ~ 1, data=data_div)
rc5k <- lm(brc.fsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.fsNA ~ fc_5km +  I(fc_5km^2), data=data_div)
rc3k <- lm(brc.fsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.fsNA ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(rc.fs.aic <- AICctab(rc0, rc5k, rc5kq, 
                      rc3k, rc3kq,
                      base=T,weights=T, logLik=T))
##       logLik AICc dLogLik dAICc df weight
## rc0    5.4   -5.5  0.0     0.0  2  0.482 
## rc3kq  8.7   -3.7  3.3     1.8  4  0.197 
## rc5k   5.9   -2.8  0.5     2.7  3  0.126 
## rc5kq  8.1   -2.5  2.7     3.0  4  0.108 
## rc3k   5.5   -2.1  0.1     3.4  3  0.087

Summary of 2nd best model, as the top is of absence of effect (~1)

summary(rc3kq)
## 
## Call:
## lm(formula = brc.fsNA ~ fc_3km + I(fc_3km^2), data = data_div)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.13558 -0.09064 -0.03247  0.07741  0.23095 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.6952262  0.2224103  -3.126   0.0122 *
## fc_3km       0.0409921  0.0173432   2.364   0.0424 *
## I(fc_3km^2) -0.0007333  0.0002926  -2.507   0.0335 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1353 on 9 degrees of freedom
## Multiple R-squared:  0.4227, Adjusted R-squared:  0.2944 
## F-statistic: 3.295 on 2 and 9 DF,  p-value: 0.08439
rc1.fs <- rc0# rename top model to save and plot below

Inspect residuals:

sim <- simulateResiduals(rc3kq)
plot(sim)

testDispersion(sim) # ok

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.82544, p-value = 0.816
## alternative hypothesis: two.sided
hnp::hnp(rc3kq)
## Gaussian model (lm object)

Prediction:

Total abundance

Modelling

### LM "gamma"  abundanec for regular model selection
gab0 <- glm.nb(ab.land_FS ~ 1, data=data_div)
gab5k <- glm.nb(ab.land_FS ~ fc_5km, data=data_div)
gab5kq <- glm.nb(ab.land_FS ~ fc_5km +  I(fc_5km^2), data=data_div)
gab3k <- glm.nb(ab.land_FS ~ fc_3km, data=data_div)
gab3kq <- glm.nb(ab.land_FS ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(gab.fs.aic <- AICctab(gab0, gab5k, gab5kq, 
                       gab3k, gab3kq,
                       base=T,weights=T, logLik=T))
##        logLik AICc  dLogLik dAICc df weight
## gab5k  -70.5  150.1   2.1     0.0 3  0.399 
## gab0   -72.6  150.5   0.0     0.5 2  0.318 
## gab3k  -71.2  151.3   1.4     1.2 3  0.217 
## gab5kq -70.4  154.5   2.2     4.5 4  0.043 
## gab3kq -71.0  155.8   1.6     5.7 4  0.023

Summary 2nd best model, as top model was the one of absence of effect (~1):

summary(gab5k)
## 
## Call:
## glm.nb(formula = ab.land_FS ~ fc_5km, data = data_div, init.theta = 4.537178159, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6496  -0.9873  -0.2873   0.6438   1.8333  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  5.93023    0.31931  18.572   <2e-16 ***
## fc_5km      -0.02485    0.01110  -2.239   0.0251 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(4.5372) family taken to be 1)
## 
##     Null deviance: 17.268  on 11  degrees of freedom
## Residual deviance: 12.411  on 10  degrees of freedom
## AIC: 147.08
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  4.54 
##           Std. Err.:  1.83 
## 
##  2 x log-likelihood:  -141.082
gab.fs <- gab5k# rename top model to save and plot below

Inspect residuals:

sim <- simulateResiduals(gab5k)
plot(sim)

testDispersion(sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.1614, p-value = 0.56
## alternative hypothesis: two.sided
hnp::hnp(gab5k)
## Negative binomial model (using MASS package)

boot::glm.diag.plots(gab5k)

# R2 based on Nakagawa et al. 2017: variação explicada pelos efeitos fixos dividida pela total
gab5k.r2 <- r2(gab5k)
gab3k.r2 <- r2(gab3k)

## build AIC table do add model selection results:
aic.tab <- rbind(aic.tab, c("Total abundance FS", 
                            "FC at 5km", 
                            round(gab.fs.aic$AICc[1],2),
                            round(gab.fs.aic$dAICc[1],1), 
                            gab.fs.aic$df[1],
                            round(gab.fs.aic$weight[1],2),
                            round(gab5k.r2$R2_Nagelkerke,2)),
                 c(" ", "~ 1 (NULL)", 
                   round(gab.fs.aic$AICc[2],2), 
                   round(gab.fs.aic$dAICc[2],1),
                   gab.fs.aic$df[2], 
                   round(gab.fs.aic$weight[2],2),
                    "-"),
                 c(" ", "FC at 3km", 
                            round(gab.fs.aic$AICc[3],2),
                            round(gab.fs.aic$dAICc[3],1), 
                            gab.fs.aic$df[3],
                            round(gab.fs.aic$weight[3],2),
                            round(gab3k.r2$R2_Nagelkerke,2)))

Prediction:

bRaup-Crick (abundance-based)

Modelling:

### LM gamma for regular model selection
rc0 <- lm(brc.abund.fsNA ~ 1, data=data_div)
rc5k <- lm(brc.abund.fsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.abund.fsNA ~ fc_5km +  I(fc_5km^2), data=data_div)
rc3k <- lm(brc.abund.fsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.abund.fsNA ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(rc.fs.aic <- AICctab(rc0, rc5k, rc5kq, 
                      rc3k, rc3kq,
                      base=T,weights=T, logLik=T))
##       logLik AICc dLogLik dAICc df weight
## rc5k   0.2    8.6  3.8     0.0  3  0.596 
## rc3k  -0.7   10.4  2.9     1.9  3  0.235 
## rc0   -3.6   12.5  0.0     3.9  2  0.085 
## rc5kq  0.3   13.1  3.9     4.5  4  0.062 
## rc3kq -0.7   15.1  2.9     6.6  4  0.022

Summary top model:

summary(rc5k)
## 
## Call:
## lm(formula = brc.abund.fsNA ~ fc_5km, data = data_div)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.39341 -0.17488 -0.01206  0.15318  0.50889 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.90123    0.17499   5.150 0.000431 ***
## fc_5km      -0.01796    0.00606  -2.964 0.014189 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2602 on 10 degrees of freedom
## Multiple R-squared:  0.4677, Adjusted R-squared:  0.4144 
## F-statistic: 8.786 on 1 and 10 DF,  p-value: 0.01419
rc2.fs <- rc5k# rename top model to save and plot below

Inspect residuals:

sim <- simulateResiduals(rc5k)
plot(sim)

testDispersion(sim) # ok

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc5k)
## Gaussian model (lm object)

Prediction

NON-FOREST SPECIALISTS

Gamma diversity

Modelling:

### LM alpha for regular model selection
g0 <- glm(gamma_NFS ~ 1, data=data_div, family= Gamma (link = "identity"))
g5k <- glm(gamma_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
g5kq <- glm(gamma_NFS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
g3k <- glm(gamma_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
g3kq <- glm(gamma_NFS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))


## AIC-based model selection
(g.nfs.aic <- AICctab(g0, 
                      g5k, g5kq, 
                      g3k, g3kq,
                      base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## g3k  -27.3   63.5   7.6     0.0 3  0.7911
## g5k  -29.2   67.5   5.7     4.0 3  0.1095
## g3kq -27.1   68.0   7.8     4.4 4  0.0864
## g5kq -29.2   72.2   5.7     8.6 4  0.0105
## g0   -34.9   75.2   0.0    11.6 2  0.0024

Summary top model:

summary(g3k)
## 
## Call:
## glm(formula = gamma_NFS ~ fc_3km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.38392  -0.22774  -0.05573   0.15121   0.43058  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 18.37185    2.13069   8.622 6.07e-06 ***
## fc_3km      -0.30280    0.05227  -5.793 0.000175 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.08050867)
## 
##     Null deviance: 2.71637  on 11  degrees of freedom
## Residual deviance: 0.78116  on 10  degrees of freedom
## AIC: 60.54
## 
## Number of Fisher Scoring iterations: 4
g.nfs <- g3k # rename top model to save and plot below

Inspect residuals:

sim <- simulateResiduals(g3k)
plot(sim)

testDispersion(sim) # a lil under

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.90839, p-value = 0.96
## alternative hypothesis: two.sided
hnp::hnp(g3k) ## for glm poisson or neg binom (maybe also others)
## Gamma model

boot::glm.diag.plots(g3k) ## for glm poisson or neg binom

Prediction:

Alpha diversity

Modelling:

### LM alpha for regular model selection
a0 <- glm(alpha_NFS ~ 1, data=data_div, family= Gamma (link = "identity"))
a5k <- glm(alpha_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
a5kq <- glm(alpha_NFS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
a3k <- glm(alpha_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
a3kq <- glm(alpha_NFS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))


## AIC-based model selection
(a.nfs.aic <- AICctab(a0, a5k, a5kq, 
                      a3k, a3kq,
                      base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## a3k  -16.7   42.3   5.7     0.0 3  0.490 
## a5k  -16.9   42.9   5.4     0.6 3  0.371 
## a5kq -16.2   46.2   6.1     3.9 4  0.071 
## a3kq -16.5   46.6   5.9     4.3 4  0.057 
## a0   -22.3   50.0   0.0     7.7 2  0.011

Summary top model:

summary(a3k)
## 
## Call:
## glm(formula = alpha_NFS ~ fc_3km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.68238  -0.15994   0.00108   0.11140   0.55352  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.82215    0.83950   6.935 4.02e-05 ***
## fc_3km      -0.09207    0.02113  -4.357  0.00143 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.1163241)
## 
##     Null deviance: 3.1391  on 11  degrees of freedom
## Residual deviance: 1.2537  on 10  degrees of freedom
## AIC: 39.336
## 
## Number of Fisher Scoring iterations: 6
a.nfs <- a3k # rename top model to save and plot below

Inspect residuals:

sim <- simulateResiduals(a3k)
plot(sim)

testDispersion(sim) # a lil under

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.93258, p-value = 0.928
## alternative hypothesis: two.sided
hnp::hnp(a3k) ## for glm poisson or neg binom (maybe also others)
## Gamma model

boot::glm.diag.plots(a3k) ## for glm poisson or neg binom

Prediction:

Beta diversity additive

Modelling

### LM Beta for regular model selection
b0 <- glm(betad_NFS~ 1, data=data_div, family= Gamma (link = "identity"))
b5k <- glm(betad_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
b5kq <- glm(betad_NFS ~ fc_5km +  I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
b3k <- glm(betad_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
b3kq <- glm(betad_NFS ~ fc_3km +  I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))

## AIC-based model selection
(b.nfs.aic <- AICctab(b0, b5k, b5kq, 
                      b3k, b3kq, 
                      base=T,weights=T, logLik=T))
##      logLik AICc  dLogLik dAICc df weight
## b3k  -25.0   58.9   5.9     0.0 3  0.753 
## b5k  -26.8   62.7   4.0     3.7 3  0.116 
## b3kq -24.6   62.9   6.3     3.9 4  0.105 
## b5kq -26.6   67.0   4.2     8.1 4  0.013 
## b0   -30.9   67.1   0.0     8.2 2  0.013

Summary top model:

summary(b3k)
## 
## Call:
## glm(formula = betad_NFS ~ fc_3km, family = Gamma(link = "identity"), 
##     data = data_div)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.71147  -0.22545  -0.03265   0.18829   0.45483  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 12.59989    1.68741   7.467 2.14e-05 ***
## fc_3km      -0.21226    0.04078  -5.205 0.000399 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.1116167)
## 
##     Null deviance: 3.2119  on 11  degrees of freedom
## Residual deviance: 1.2312  on 10  degrees of freedom
## AIC: 55.933
## 
## Number of Fisher Scoring iterations: 4
b.nfs <- b3k # rename top model to save and plot below

Inspect residuals:

sim <- simulateResiduals(b3k)
plot(sim)

testDispersion(sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.65877, p-value = 0.68
## alternative hypothesis: two.sided
hnp::hnp(b3k) ## for glm poisson or neg binom (maybe also others)
## Gamma model

boot::glm.diag.plots(b3k) ## for glm poisson or neg binom

Prediction:

bRaup-Crick (presence/absence-based)

Modelling

### LM gamma for regular model selection
rc0 <- lm(brc.nfsNA ~ 1, data=data_div)
rc5k <- lm(brc.nfsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.nfsNA ~ fc_5km +  I(fc_5km^2), data=data_div)
rc3k <- lm(brc.nfsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.nfsNA ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(rc.nfs.aic <- AICctab(rc0, rc5k, rc5kq, 
                       rc3k, rc3kq,
                       base=T,weights=T, logLik=T))
##       logLik AICc dLogLik dAICc df weight
## rc3k  -0.5   10.0  2.3     0.0  3  0.393 
## rc0   -2.8   11.0  0.0     1.0  2  0.244 
## rc5k  -1.1   11.2  1.7     1.2  3  0.220 
## rc3kq  0.6   12.5  3.4     2.4  4  0.116 
## rc5kq -0.8   15.4  2.0     5.3  4  0.027

Summary of 2nd best model, 1st wqs of absence of effect (~1):

summary(rc3k)
## 
## Call:
## lm(formula = brc.nfsNA ~ fc_3km, data = data_div)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32362 -0.21530 -0.05601  0.23870  0.39743 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.32268    0.18990   1.699   0.1201  
## fc_3km      -0.01402    0.00647  -2.167   0.0555 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2766 on 10 degrees of freedom
## Multiple R-squared:  0.3195, Adjusted R-squared:  0.2515 
## F-statistic: 4.696 on 1 and 10 DF,  p-value: 0.05545
rc1.nfs <- rc3k # rename top model to save and plot below

Inspect residuals:

sim <- simulateResiduals(rc3k)
plot(sim)

testDispersion(sim) # a lil under

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc3k)
## Gaussian model (lm object)

Prediction:

Total abundance

Modelling:

### LM "gamma"  abundanec for regular model selection
gab0 <- glm.nb(ab.land_NFS ~ 1, data=data_div)
gab5k <- glm.nb(ab.land_NFS ~ fc_5km, data=data_div)
gab5kq <- glm.nb(ab.land_NFS ~ fc_5km +  I(fc_5km^2), data=data_div)
gab3k <- glm.nb(ab.land_NFS ~ fc_3km, data=data_div)
gab3kq <- glm.nb(ab.land_NFS ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(gab.nfs.aic <- AICctab(gab0, gab5k, gab5kq, 
                        gab3k, gab3kq,
                        base=T,weights=T, logLik=T))
##        logLik AICc  dLogLik dAICc df weight
## gab5k  -63.9  136.9   5.3     0.0 3  0.404 
## gab3k  -64.3  137.6   4.9     0.7 3  0.284 
## gab5kq -62.1  137.8   7.2     0.9 4  0.252 
## gab3kq -63.7  141.2   5.5     4.3 4  0.047 
## gab0   -69.2  143.8   0.0     6.9 2  0.013

Summary top model:

summary(gab5k)
## 
## Call:
## glm.nb(formula = ab.land_NFS ~ fc_5km, data = data_div, init.theta = 4.164123767, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.8516  -0.4183   0.2175   0.4846   0.8634  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  5.93704    0.33661  17.638  < 2e-16 ***
## fc_5km      -0.04838    0.01182  -4.094 4.24e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(4.1641) family taken to be 1)
## 
##     Null deviance: 29.266  on 11  degrees of freedom
## Residual deviance: 12.765  on 10  degrees of freedom
## AIC: 133.88
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  4.16 
##           Std. Err.:  1.77 
## 
##  2 x log-likelihood:  -127.88
gab.nfs <- gab5k # rename top model to save and plot below

Inspect residuals:

sim <- simulateResiduals(gab5k)
plot(sim)

testDispersion(sim) # a lil under

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.45825, p-value = 0.288
## alternative hypothesis: two.sided
hnp::hnp(gab5k)
## Negative binomial model (using MASS package)

boot::glm.diag.plots(gab5k)

# R2 based on Nakagawa et al. 2017: variação explicada pelos efeitos fixos dividida pela total
gab5k.r2 <- r2(gab5k)
gab5kq.r2 <- r2(gab5kq)
gab3k.r2 <- r2(gab3k)

## build AIC table do add model selection results:
aic.tab <- rbind(aic.tab, c("Total abundance NFS", 
                            "FC at 5km", 
                            round(gab.nfs.aic$AICc[1],2),
                            round(gab.nfs.aic$dAICc[1],1), 
                            gab.nfs.aic$df[1],
                            round(gab.nfs.aic$weight[1],2),
                            round(gab5k.r2$R2_Nagelkerke,2)),
                 c(" ", "FC at 3km",
                   round(gab.nfs.aic$AICc[2],2), 
                   round(gab.nfs.aic$dAICc[2],1),
                   gab.nfs.aic$df[2], 
                   round(gab.nfs.aic$weight[2],2),
                   round(gab3k.r2$R2_Nagelkerke,2)),
                 c(" ", "Non-linear FC at 5km",
                   round(gab.nfs.aic$AICc[3],2), 
                   round(gab.nfs.aic$dAICc[3],1),
                   gab.nfs.aic$df[3], 
                   round(gab.nfs.aic$weight[3],2),
                   round(gab5kq.r2$R2_Nagelkerke,2)
                 ))

Prediction:

bRaup-Crick (abundance-based)

Modelling

### LM gamma for regular model selection
rc0 <- lm(brc.abund.nfsNA ~ 1, data=data_div)
rc5k <- lm(brc.abund.nfsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.abund.nfsNA ~ fc_5km +  I(fc_5km^2), data=data_div)
rc3k <- lm(brc.abund.nfsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.abund.nfsNA ~ fc_3km +  I(fc_3km^2), data=data_div)

## AIC-based model selection
(rc.nfs.aic <- AICctab(rc0, rc5k, rc5kq, 
                       rc3k, rc3kq,
                       base=T,weights=T, logLik=T))
##       logLik AICc dLogLik dAICc df weight
## rc3k  -7.5   24.0  1.9     0.0  3  0.377 
## rc0   -9.4   24.2  0.0     0.2  2  0.341 
## rc5k  -8.0   25.1  1.4     1.1  3  0.219 
## rc3kq -7.4   28.5  2.1     4.5  4  0.041 
## rc5kq -8.0   29.7  1.5     5.7  4  0.022

Summary 2nd best model, 1st is of absence of effect (~1):

summary(rc3k)
## 
## Call:
## lm(formula = brc.abund.nfsNA ~ fc_3km, data = data_div)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.25085 -0.18926  0.04247  0.27112  0.52852 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.84365    0.33999   2.481   0.0325 *
## fc_3km      -0.02260    0.01158  -1.951   0.0797 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4951 on 10 degrees of freedom
## Multiple R-squared:  0.2756, Adjusted R-squared:  0.2032 
## F-statistic: 3.805 on 1 and 10 DF,  p-value: 0.07965
rc2.nfs <- rc3k # rename top model to save and plot below

Inspect residuals:

sim <- simulateResiduals(rc3k)
plot(sim)

testDispersion(sim) # a lil under

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc3k)
## Gaussian model (lm object)

Prediction:

AICc-model selection table

Table S2. AICc-based selected models table containing the models of gamma (γ), alpha (ɑ), additive beta (β) diversities followed by abundance at landscape total, abundance at mean site and β-diversity calculated from the Raup-Crick null model approach, for both presence/absence (βRC) and abundance-based (βRC-abundance) that scored below ΔAICc=2. All tested as a function of Atlantic Forest loss for each group of species, first models separated into forest specialists species (FS; Table 1a), and followed by non-forest species (NFS; Table 1b). We show selected models ordered from the lowest to the highest AICc values ΔAICc≤ 2. We present information regarding the diversity metric tested (Diversity), the model terms for forest loss predictors (Model), the Akaike information criterion corrected for small samples (AICc), the AICc difference from the first-ranked model (ΔAICc), degrees of freedom (df) Akaike weight (ωi) and model variance explained (r²).
Diversity Model AICc dAICc df weight r2
Gamma FS FC at 5km 64.84 0 3 0.73 0.55
Alpha FS FC at 5km 47.61 0 3 0.52 0.38
FC at 3km 49.54 1.9 3 0.2 0.27
~ 1 (NULL) 49.58 2 2 0.19
Beta FS FC at 5km 53.26 0 3 0.69 0.58
RCbeta (P/A-based) FS ~ 1 (NULL) -5.48 0 2 0.48
Non-linear FC at 3km -3.7 1.8 4 0.2 0.29
Total abundance FS FC at 5km 150.08 0 3 0.4 0.44
~ 1 (NULL) 150.54 0.5 2 0.32
FC at 3km 151.31 1.2 3 0.22 0.32
RCbeta (abund-based) FS FC at 5km 8.56 0 3 0.6 0.41
FC at 3km 10.41 1.9 3 0.24 0.32
Gamma NFS FC at 3km 63.54 0 3 0.79 0.74
Alpha NFS FC at 3km 42.34 0 3 0.49 0.63
FC at 5km 42.89 0.6 3 0.37 0.61
Beta NFS FC at 3km 58.93 0 3 0.75 0.65
RCbeta (P/A-based) NFS FC at 3km 10.02 0 3 0.39 0.25
~ 1 (NULL) 10.97 1 2 0.24
FC at 5km 11.18 1.2 3 0.22 0.18
Total abundance NFS FC at 5km 136.88 0 3 0.4 0.82
FC at 3km 137.58 0.7 3 0.28 0.79
Non-linear FC at 5km 137.82 0.9 4 0.25 0.93
RCbeta (abund-based) NFS FC at 3km 24 0 3 0.38 0.2
~ 1 (NULL) 24.2 0.2 2 0.34
FC at 5km 25.09 1.09 3 0.22 0.13

FIGURES

Gamma hab

### Gamma ####
# r2(g.fs) # 0.55
## prediction g2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(g.fs, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

# r2(g.nfs) # 0.73

## prediction g1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds2 <- predict(g.nfs, type="response", newdata=new2, se.fit=T)

new2$pred <- preds2$fit
new2$ic.up <- preds2$fit + preds2$se.fit*1.96 
new2$ic.low <- preds2$fit - preds2$se.fit*1.96 

## Plot gamma richness ~ fc %
(g.hab <- ggplot(data= new, aes(x= fc_5km, y = pred)) + 
    geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
    geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
    geom_point(data= data_div, aes(x=fc_5km, y=gamma_FS), col= "#0072B2") + 
    geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
    geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D55E00") +
    geom_point(data= data_div, aes(x=fc_3km, y=gamma_NFS), col= "#D55E00") +
    xlab(" ") + ylab("Gamma") + 
    theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') + 
    ylim(0,max(data_div[2:7])+0.5) + 
    # scale_x_reverse() +
    annotate("text", y=20.5, x=38,colour= "#0072B2",
             label=expression(paste(~ R^2 ~ "= 0.55"))) +
    annotate("text", y=19, x=40, colour= "#D55E00",
             label=expression(paste(~ R^2 ~ "= 0.73"))) + labs(tag="a"))

Alpha hab

# r2(a.fs) # 0.38

## prediction a2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(a.fs, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

# r2(a.nfs) # 0.63
## prediction a1
new2 <-  data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(a.nfs, type="response", newdata=new2, se.fit=T)

new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96 
new2$ic.low <- preds$fit - preds$se.fit*1.96 

## Plot Alpha richness ~ fc_3km
(a.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) + 
  geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2", 
            alpha= .5, linetype= "dashed") +
  # geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
  geom_point(data=data_div, aes(x=fc_5km, y=alpha_FS), col= "#0072B2") +
  geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D55E00") +
  geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
  geom_point(data=data_div, aes(x=fc_3km, y=alpha_NFS), col= "#D55E00") +
  xlab("") + ylab("Alpha") + 
  scale_color_manual(breaks=c("FS", "NFS"),
                     values=c('FS'='#0072B2', 'NFS'='#D55E00')) +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=14,face="bold"),
        legend.position = 'bottom') + 
  ylim(0,max(data_div[2:7])+0.5) +
    # scale_x_reverse() +
  annotate("text", y=13, x=38,colour= "#0072B2",
           label=expression(paste(~ R^2 ~ "= 0.38")))+
  annotate("text", y=11, x=40, colour= "#D55E00",
           label=expression(paste(~ R^2 ~ "= 0.63"))) +
  labs(tag="b"))

Beta hab

### Beta add  ####
# r2(b.fs) # 0.58

## prediction b1
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(b.fs, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

# r2(b.nfs) # 0.65

## prediction b1
new2 <-  data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(b.nfs, type="response", newdata=new2, se.fit=T)

new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96 
new2$ic.low <- preds$fit - preds$se.fit*1.96 

## Plot Beta Diversity ~ fc %
(b.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) + 
  geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
  geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
  geom_point(data=data_div, aes(x=fc_5km, y=betad_FS), col= "#0072B2") +
  geom_line(data= new2, aes(x= fc_3km, y= pred), 
            col= "#D55E00") +
  geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
  geom_point(data=data_div, aes(x=fc_3km, y=betad_NFS), col= "#D55E00") +
  xlab(" ") + ylab("Beta") + 
  theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') + 
  ylim(0,max(data_div[2:7])+0.5) + 
    # scale_x_reverse() +
  annotate("text", y=18, x=40,colour= "#0072B2",
           label=expression(paste(~ R^2 ~ "= 0.58")))+
  annotate("text", y=16, x=38, colour= "#D55E00",
           label=expression(paste(~ R^2 ~ "= 0.65"))) +
  labs(tag="c"))

g.hab + a.hab + b.hab

ggsave("figures/plot.hab.jpeg",units="cm", width=20, height=8, dpi=300)

Total abundance FS & NFS

## FS
# r2(gab.fs) # 0.44

## prediction gab1
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(gab.fs, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

## NFS
# r2(gab.nfs) # 0.82

## prediction gab2
new2 <- new
preds <- predict(gab.nfs, type="response", newdata=new2, se.fit=T)

new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96 
new2$ic.low <- preds$fit - preds$se.fit*1.96 

## Plot gamma richness ~ forest cover 5km
(gab.hab <- ggplot(new, aes(x=fc_5km, y=pred)) + 
  geom_line(linetype= "dashed", col= "#0072B2", aes(alpha= 0.1)) +
  # geom_ribbon(aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1,  fill= "#0072B2") +
  geom_point(data=data_div, aes(x=fc_5km, y=ab.land_FS), col= "#0072B2") +
  geom_line(data= new2, aes(x=fc_5km, y=pred), col= "#D55E00") +
  geom_ribbon(data= new2, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.2, fill= "#D55E00") +
  geom_point(data= data_div, aes(x=fc_5km, y=ab.land_NFS), col= "#D55E00") +
    scale_x_continuous(limits= c(10,50)) +
  xlab(" ") + ylab("Total abundance") + 
  theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') + 
  # scale_x_reverse() +
    annotate("text", y=300, x=38, alpha= .5, colour= "#0072B2",
             label=expression(paste(~ R^2 ~ "= 0.44"))) +
    annotate("text", y=270, x=40, colour= "#D55E00",
           label=expression(paste(~ R^2 ~ "= 0.82"))) +
  labs(tag="d"))

### Beta RC FS & NFS

## FS
# ## NULL ## r2(rc1.fs) 

## prediction rc2.2 and rc2
new <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=30))
preds <- predict(rc1.fs, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

## NFS
# r2(rc1.nfs) # 0.25

## prediction rc2
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(rc1.nfs, type="response", newdata=new2, se.fit=T)

new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96 
new2$ic.low <- preds$fit - preds$se.fit*1.96 

## Plot RCbeta abundance ~ forest cover at 3
(rc.hab <- ggplot(data= new, aes(x=fc_3km, y=pred)) + 
  geom_line(data= new,aes(x=fc_3km, y=pred), 
            col= "#0072B2", linetype= "dashed", alpha= .2) +
  # geom_ribbon(data= new, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") + 
  geom_point(data=data_div, aes(x=fc_3km, y=brc.fsNA), col= "#0072B2") +
  geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D33E00", linetype= "dashed", alpha= .2) +
  # geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D33E00") +
  geom_point(data=data_div, aes(x=fc_3km, y=brc.nfsNA), col= "#D33E00") +
  xlab("Percent forest cover") + ylab(c(expression(bold("β"["RC"])))) + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14,face="bold"),
        legend.position = 'none') + ylim(-1,1) + 
    # scale_x_reverse() +
  annotate("text", y=0.9, x=43, alpha= .5, colour= "#D33E00",
           label=expression(paste(~ R^2 ~ "= 0.25"))) +
  # annotate("text", y=0.7, x=43, alpha= .5, colour= "#0072B2",
           # label=expression(paste(~ R^2 ~ "= 0.30"))) +
  # annotate("text", y=0.9, x=43, colour= "#D33E00",
  #          label=expression(paste(beta, "(SE) = -0.06(0.03); " ~ R^2 ~ "= 0.37"))) +
  # annotate("text", y=0.8, x=43, colour= "#0072B2",
  #          label=expression(paste(beta, "(SE) = -0.06(0.03); " ~ R^2 ~ "= 0.33"))) +
  labs(tag="e"))

Beta RC abund FS & NFS

# r2(rc2.fs) # 0.41

## prediction rc2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(rc2.fs, type="response", newdata=new, se.fit=T)

new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96 
new$ic.low <- preds$fit - preds$se.fit*1.96 

## NFS
# r2(rc2.nfs) # 0.21

## prediction rc1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(rc2.nfs, type="response", newdata=new2, se.fit=T)

new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96 
new2$ic.low <- preds$fit - preds$se.fit*1.96 

## RCbeta abund ~ FC%
(rc.ab.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) + 
  geom_line(data= new,aes(x=fc_5km, y=pred), col= "#0072B2") +
  geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") + 
  geom_point(data=data_div, aes(x=fc_5km, y=brc.abund.fsNA), col= "#0072B2") +
  geom_line(data= new2, aes(x=fc_3km, y=pred), linetype= "dashed", col= "#D55E00", alpha=.2) +
  # geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
  geom_point(data=data_div, aes(x=fc_3km, y=brc.abund.nfsNA), col= "#D55E00") +
  xlab(" ") + ylab(c(expression(bold("β"["RC-abundance"])))) + 
  theme(axis.text=element_text(size=14), axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
    ylim(-1,1.1) + 
    # scale_x_reverse() +
  annotate("text", y=0.9, x=45, colour= "#0072B2",
           label=expression(paste(~ R^2 ~ "= 0.41"))) +
  annotate("text", y=0.75, x=43, alpha= .5, colour= "#D55E00",
           label=expression(paste( ~ R^2 ~ "= 0.21"))) +
  labs(tag="f"))

Final plot

(plot.hab <- (g.hab + a.hab + b.hab) /
    (gab.hab + rc.hab + rc.ab.hab))

ggsave("figures/figure2.jpeg",units="cm", width=28, height=15, dpi=300)